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Abstract It has been proposed that the milhon degree temperature of the corona 
is due to the combined effect of barely-detectable energy releases, so called 
nanoflares, that occur throughout the solar atmosphere. Alas, the nanoflare 
density and brightness implied by this hypothesis means that conclusive verifi- 
cation is beyond present observational abilities. Nevertheless, we investigate the 
plausibility of the nanoflare hypothesis by constructing a magnetohydrodynamic 
(MHD) model that can derive the energy of a nanoflare from the nature of an 
ideal kink instability. The set of energy-releasing instabilities is captured by 
an instability threshold for linear kink modes. Each point on the threshold is 
associated with a unique energy release and so we can predict a distribution of 
nanoflare energies. When the linear instability threshold is crossed, the instability 
enters a nonlinear phase as it is driven by current sheet reconnection. As the 
ensuing flare erupts and declines, the field transitions to a lower energy state, 
which is modelled by relaxation theory, i.e., helicity is conserved and the ratio 
of current to field becomes invariant within the loop. We apply the model so 
that all the loops within an ensemble achieve instability followed by energy- 
releasing relaxation. The result is a nanoflare energy distribution. Furthermore, 
we produce different distributions by varying the loop aspect ratio, the nature of 
the path to instability taken by each loop and also the level of radial expansion 
that may accompany loop relaxation. The heating rate obtained is just sufficient 
for coronal heating. In addition, we also show that kink instability cannot be as- 
sociated with a critical magnetic twist value for every point along the instability 
threshold. 
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1. Introduction 

The theory of nanoflare coronal heating (Parker, 1988) postulates that small 
flares are sufficiently numerous to maintain a coronal temperature of millions 
of degrees. Observational studies have so far been unable to show conclusively 
whether nanoflares occur frequently enough to heat the corona (Krucker and 
Benz, 1998; Parnell and Jupp, 2000; Aschwanden and Parnell, 2002; Parnell, 
2004): the smaller a flare the harder it is to distinguish from the coronal back- 
ground. Nevertheless, we propose a coronal loop model, the purpose of which is 
to test the viability of the nanoflare theory with respect to coronal heating. 

Random convective motions at the photosphere increase the energy of the 
magnetic fields that define coronal loops (plasma beta, /3«0.01). The magnetic 
field is repeatedly twisted by such motions; coronal conductivities arc so large 
that the magnetic flux and plasma can be regarded as being frozen together. 
We assume that the kinetic energy imparted by the photosphere is dissipated 
via Direct Current heating: the timescale for photospheric turbulence is long 
compared to the Alfven time (Klimchuk, 2006). Hence, the loop can be repre- 
sented as moving through a series of force-free states: V x B = a{f)B, where 
a = {noj ■ B)/{\B\'^) is the ratio of current density to magnetic field and r is a 
position vector (Woltjer, 1958). Magnetic stresses build within the loop until an 
instability is reached or dissipation is otherwise triggered. 

Magnetic reconnection is thought to be the mechanism by which the excess 
magnetic energy is converted into heat. Observations of large-scale flares have 
revealed stong evidence for such a process (Fletcher, 2009; Qiu, 2009). In ad- 
dition, several 3D MHD models have shown that current sheets form during 
the nonlinear phase of an ideal kink instability (Baty and Heyvaerts, 1996; 
Velli, Lionello, and Einaudi, 1997; Arber, Longbottom, and Van der Linden, 
1999; Baty, 2000). The kink instability gives rise to helical current sheets that 
become the site of Ohmic dissipation, i.e., a heating event. These models are 
restricted to a narrow range of initial loop conflgurations. Expanding this range 
so that one could determine the relationship between say, the initial a-proflle and 
the resulting energy release would be too computationally expensive. However, 
the energy release can be calculated without recourse to following the complex 
dynamics of the reconnection process (Heyvaerts and Priest, 1984). 

When a magnetic fleld becomes unstable, relaxation theory predicts the field 
will transition to the lowest energy state that conserves total magnetic axial 
flux and global magnetic helicity (Taylor, 1974, 1986). The minimum energy (or 
relaxed) state is the well-known constant a or linear force-free field, \7xB ^ aB. 
The helicity (K) measures the degree to which the magnetic field is linked with 
itself (Berger, 1999). However, the relative helicity (Berger and Field, 1984; Finn 
and Antonsen, 1985) is more useful since it is gauge invariant: 

K= [ {A + A')-{B-B') dV, (1) 
Jv 

where A is the magnetic potential, B' is the potential field with the same 
boundary conditions and A' is the corresponding vector potential. 
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Helicity conservation is not absolute. During relaxation, lielicity is still subject 
to global resistive diffusion, but the change in helicity is negligible when com- 
pared to the drop in magnetic energy, so long as dissipation predominantly occurs 
within thin current sheets. The rates of dissipation for helicity and magnetic 
energy {W) are 



where j=V x i?//io is the current density, I is the length scale of magnetic 
variation (i.e., current sheet thickness), L is the global length scale and rj is the 
resistivity (Browning, 1988). Using K « B^L and W ~ B"^ /2^q, the ratio 
of the dissipation rates reduces to l/L. Hence, dtK/K <^ dtW/W ii I ^ L, 
which is expected to be well satisfied for reconnecting current sheets within the 
highly conductive corona, where global resistive diffusion of helicity and energy 
are negligible. The relative sizes of the dissipation rates have been confirmed 
by MHD simulations, despite the coarseness of numerical grids (the difference 
between dissipation rates becomes more pronounced as the resistivity becomes 
smaller and falls below numerical precision). Browning et al. (2008) showed that 
during the relaxation of a marginally (kink) unstable loop SK/K ~ 10"'' and 
S W/W ^10^^. Detailed estimates of coronal helicity dissipation are given by 
Berger (1984); further justification for helicity-conserving relaxation is provided 
by laboratory experiments (Hcidbrink and Dang, 2000; Taylor, 1986) 

Helicity conservation combined with the invariant nature of the relaxed a- 
profile imply that helicity has simply become more evenly distributed within 
the loop. Thus, the relaxed a can be calculated, as can the amount of energy 
liberated from the magnetic field during relaxation. The latter quantity can 
be interpreted as an upper limit to the heating event energy, since complete 
relaxation may not be attained. Issues concerning plasma response are outside 
the scope of the model presented in this paper. 

Repeated episodes of slow photospheric driving may trigger an ideal MHD in- 
stability. Ideal (not resistive) instabilities are required in order for the time scale 
of the instability to match observations of the highly conducting corona, where 
resistive time scales are very long. Browning and Van der Linden (2003) describe 
how a dynamic heating event is caused whenever the field exceeds the threshold 
for a linear kink instability in a cylinderical coronal loop. Extensive numerical 
simulations (Galsgaard and Nordlund, 1997; Velli, Lionello, and Einaudi, 1997; 
Lionello et al., 1998; Baty, 2000; Gerrard et al., 2001) have demonstrated how 
energy release occurs during the nonlinear phase of the instability. Relaxation 
of the field towards a constant-a state has also been demonstrated (Browning et 
al., 2008; Hood, Browning, and Van der Linden, 2009) alongwith a measurement 
of the energy release that is in good agreement with relaxation theory. 

In summary, the model presented here is based on the idea that coronal 
loops are moved by photospheric motions through a series of force-free equilibria 





(3) 
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that will eventually culminate in a heating event whenever the ideal instability 
threshold is crossed. Bareford, Browning and Van der Linden (2010, hereafter 
BBV2010) calculated the heating event distributions for an ensemble of loops 
that possessed net current. This was done using a simple cylinderical field model 
in which the current profile {a{r)) of the stressed field is represented by a two pa- 
rameter family. (Theoretically, only two states on the closed instability threshold 
for two-parameter loops could yield zero- net-current configurations.) The result- 
ing distributions were compatible with the energies required for coronal heating. 
Following BBV2010, we use a simple cylindrical field model, and represent the 
nonlinear force-free field using a three-parameter family of current profiles in 
which a is a piecewisc constant. Wc now attempt to improve the realism of the 
model by including a mechanism that requires each loop to have zero net current. 
Thus, we consider the effects of twisting motions localised within the loop cross 
section (Hood, Browning, and Van der Linden, 2009). Previously, the random 
nature of photospheric motions was represented by allowing different parts of 
the loop interior to vary independently. It is more reasonable however to assume 
some level of correlation, since it is likely that the whole of a loop footpoint will 
be subjected to the same convective eddy. Furthermore, we also consider the 
consequence of varying the aspect ratio of the loop [L / R). This paper will also 
check to see if there is a twist-based parameter that can be used as a simple 
diagnostic for loop instability, e.g.. Hood and Priest (1979). 

The composition of the loop model and the equations used to express the 
magnetic field are given in the next section. Section 3 describes how the loop 
instability threshold is calculated and presents the equations used to determine 
the energy release associated with each relaxation. This section also analyses the 
loop configurations that follow the instability threshold. In addition, the loop's 
path to instability is explained. The nanofiarc energy distributions produced by 
the model arc presented in Section 4. Finally, in the last section, the results are 
summarised and discussed. 



2. Model 

A loop is considered to evolve through equilibria as it is driven by slow pho- 
tospheric footpoint motions. An idealised model of a straight cylindrical loop 
is used with the photosphere represented by two planes at z = 0, L; however, 
the essential physics should apply to more complex geometries. The stressed 
field is line-tied with (in general) a non-uniform a{r) (where a = ^oj^^/B). 
This is represented by a three-parameter family of piecewise-constant-a profiles. 
The model discussed here (Figure 1) is an improved version of the one used by 
BBV2010, which allowed loops to have net current (i.e., an azimuthal field was 
present in the potential envelope). A net current contradicts the idea that the 
magnetic field is twisted by convective motions local to the loop footpoints. If 
the twisting motions are confined to some localised region, the field surrounding 
the loop should be purely axial. The currents generated by the twisting of the 
fields within the loop should close locally such that the loop carries zero net 
current. Hence, a current neutralisation layer is introduced here, defined such 
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Figure 1. Schematic of a straightened coronal loop in the r-6 plane (left) and in the 
r-z plane (right). The loop, comprises a core (dark grey), an outer layer (hght grey) 
and a current neutralisation layer (blue); the whole loop is embedded in a potential 
envelope (white). The core radius is half the loop radius and 1/6 the envelope radius 
{Ri:R2.Ri:RA = 0.5:0.9:1:3). The loop aspect ratio (i / R-i) in this figure is 20. 

that the azimuthal field {Bg) always falls to zero at the loop boundary {R-i). 
Note that Hood, Browning and Van der Linden (2009) undertook 3D numerical 
simulations with initial fields taken as twisted states with zero net current. There 
are, however, only two such fields on the marginal instability curve of BBV2010, 
for which the current due to a2 cancels that due to ai (that is ai ~ =f2.48, 
a2 ~±0.95). Here, we construct a whole family of current-neutralised fields by 
adding the extra layer. 

Following on from BBV2010 and Browning and Van der Linden (2003), it 
is proposed that a relaxation event is triggered when the loop's field becomes 
linearly unstable. The energy released, due to fast magnetic reconnection dur- 
ing the nonlinear development of the instability, can then be calculated using 
relaxation theory. 

The loop's radial a-profilc is approximated by a piecewise-constant function 
featuring three parameters. The ratio of current to magnetic field is ai in the 
core, a2 in the outer layer, in the neutralisation layer and zero in the potential 
envelope. The free parameters are ai and 0^2, whereas is dependent on the 
first two and is determined by the requirement of zero net current. Note that the 
magnetic field is continuous everywhere, whereas the current has discontinuities. 
Recent work indicates that a discontinuous current has little discernable effect 
on the dynamics when compared to similar but continuous a-profiles (Hood, 
Browning, and Van der Linden, 2009). Following the investigations of Brown- 
ing et al. (2008), the outer surface of the potential envelope, representing the 
background corona, is placed at i?4 = 3 (thrice the loop radius). 

The fields are expressed in terms of the well-known Bessel function model, 
generalised to the concentric layer geometry (Melrose, NichoUs, and Brodcrick, 
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1994; Browning and Van der Linden, 2003; Browning et al., 2008; BBV2010). 
Thus, the field equations for the four regions (core, outer layer, neutralisation 
layer and envelope) are as follows: 

Si, - SiJo(|«i|r) (4) 
Bie = aiBiM\ai\r) < r < Ri (5) 

B2. = B2M\a2\r) + C2Yoi\a2\r) (6) 
B2e = <J2{B2Ji{\a2\r) + C2Yi{\a2\r)) Ri < r < R2 (7) 

Biz = BsMla^lr) + CsYoilaslr) (8) 
Bse ^ ^3{B3Ji{\a3\r) + C3Yi{\a3\r)) R2 < r < R3 (9) 

Biz = Bi (10) 
B4e = i?3<r<i?4 (11) 

where ai = (z = 1,2,3) represent the sign of a^. The fields must be continuous 
at the inner radial boundaries, Ri, R2 and R3. (We choose i?i=0.5, i?2 = 0.9 
and i?3 = 1, so that most of the loop is similar to previous work, BBV2010, with 
a thin current neutralisation layer between R2 and R3) Therefore, the constants 
Bj and Cj {j = 2,3,4) can be expressed like so: 



B. = 



n\a2\B iR 
2 

7T\a2\BiR 



-('Ji,2Jii\ai\Ri)Yoi\a2\Ri) - Joi\ai\Ri)Yii\a2\Ri)) (12) 



C2 = 1^1 (^j^(^\a,\R,)J,{\a2\Ri) - ^1,2 Ji (|«i |fli) Jo(|a2 |i?i)) (13) 

B3 = ^^^^^^(fT2,3i^i(Ia2|i?2)>o(|a3|i?2)-Fo(|a2|i?2)yi(|«3|i?2))(14) 
C3 = ^^^4^f^^o(|a2|i?2)Ji(|a3|i?2)-^2,3Fi(|a2|i?2)Jo(|a3|i?2)) (15) 



2 

Bi = BsGoila^lRs) (16) 

C4 = (17) 
where 

Fo,ii^) = Jo.iix) + ^Yo.iix) (18) 
^2 

Go.i(x) = Jo,i(x) + ^ro.i(x) (19) 

-D3 
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The value of 0:3 (the ncutraUsation layer current) is found, for a given (ai, a2), 
by numerical solution of i?30(i?3) = 0, ensuring that the net current is zero and 
that the azimuthal field vanishes outside the loop, see Equation 9. Thus, the 
equilibrium parameter space remains 2D (i.e., it is determined by ai, Q!2) - the 
field profiles for a selection of loop configurations are given in Appendix B. 
The magnetic flux through the loop and envelope is conserved: 



/ 2Trr*B;dr* = —-^M\ai\Ri) 

Jo m\ 



+ ^^^^^^i(|«2|i?2) - 2^i?iJi(|ai|i?i) ( fi^ ) 

+ ^^(i?3Gi(|a3|i?3)-i?2Gi(|a3|i?2)) 

+ TrBi(^Rl-Ri^ (20) 

where the asterisks denote dimensionless quantities. Hence, in the model, ip* is 
normalised to 1 and Bi is determined (noting that, in Equations 12-17, Bj and 
Cj are functions of Bi). The normalised coronal loop radius (i?3 = l) is itself 
used to normalise the loop length (e.g. L = 20iZ3), see Figure 1. 

We assume that the loop evolves through a sequence of two- alpha fields (Equa- 
tions 4-11) as it is twisted by photospheric footpoint motions. It has been shown 
that the introduction of magnetic twist gives the coronal loop a circular cross 
section (Klimchuk, Antiochos, and Norton, 2000). The loop model presented here 
has the same cross-sectional shape, but the loop radius {R3) is held constant 
throughout the simulated photospheric driving. Purely azimuthal photospheric 
motions would cause a small expansion of the loop (Browning and Hood, 1989) 
which we ignore; alternatively, small radial footpoint motions must be allowed 
in order to maintain constant loop radius. In any case, the sequence of loop 
equilibria explored by our model is clearly a small subset of the possible variation 
in field profiles that might arise from photospheric motions. As these random 
motions proceed, the loop continues to evolve through force-free equilibria until 
it becomes linearly unstable. We now discuss the calculation of the instability 
onset. 



3. Linear kink instability thresholds 

A coronal loop's instability is constrained by the line-tying of the photospheric 
footpoints (Hood, 1992). Hence, all perturbations are required to vanish at the 
loop ends (2 = 0, L). Any linear perturbation can be decomposed as a sum, 

^ /(r, z)e*"'^e'''*. We need only consider the m~l term however, since this 

m=0 

azimuthal mode has been found to be the least stable (Van der Linden and 
Hood, 1999) and is the dominant instability. The effect of such perturbations 
on the coronal loop are represented by the standard set of linearised ideal MHD 
equations. When the growth rate of a perturbation transitions from a negative 
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value to a positive one, the loop has reached the threshold of an ideal linear 
instability. The instability threshold is a curve in 2-dimensional a-space (ai, 
02). The properties of the loop (e.g., ai, 0,2 and a^) at these threshold points 
can be found by substituting the perturbation function into the linearised MHD 
equations, leading to an eigenvalue equation for the growth rates (Priest, 1987). 
The growth rates and cigcnfunctions of the most unstable modes are found 
numerically, for line-tied fields, with the CILTS code, described in Browning and 
Van dor Linden (2003) and Browning ct al. (2008). CILTS can be configured such 
that one of the loop's free a parameters is fixed whilst the other is incremented. 
The code terminates as soon as the real part of the cigenfunction falls below 
zero, i.e., the loop is no longer unstable to kink perturbations. 




unstable 



-6-4-2 2 4 6 

a, 

Figure 2. The instability thresholds for L / Ra = 10 (red), 20 (black) and 30 (blue). A 
closed stability region is formed by the reversal lines (dashed). 

Figure 2 shows the instability threshold curves mapped by the CILTS code 
for three values of loop aspect ratio (L /i?3 ~ 10, 20, 30). The longer the loop the 
smaller the a-value required for instability, since, if the radius is held constant, 
longer loops are less affected by line-tying stabilisations (Priest and Hood, 1979). 
The addition of a current neutralisation layer prevents the threshold curves from 
closing near the a-space origin; this is unlike the net current case for which the 
threshold is a closed curve (BBV2010, Figure 2). The open shape is indeed similar 
to the instability curve for loops with a conducting wall at (Browning and Van 
der Linden, 2003). This is because the eigenfunctions almost vanish at i?3, see 
Figures 4-6. Also, if ai is small, will be opposite in sign to a2, and the outer 
layer is stabilised by the neutralisation layer. However, the loop configurations 
become unrealistic as we increase the magnitude of a2 and the axial field reverses. 
Positions outside the reversal lines (Figure 2) represent loops that have axial 
fields of mixed polarity. These configurations cannot represent states attained 
by the twisting of a unipolar loop. 
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20 40 60 80 100 

"1 

Threshold Point 



Figure 3. On the left is the right half (i.e., where qi > 0) of the instability threshold for 
a loop of aspect ratio 20. On the right is the variation in along the ID representation 
of the instability threshold. Notice that the tic marks along the Threshold Point axis 
correspond with the numbers that follow the threshold curves shown in the left plot. 

Before proceeding to ealculate the energy release properties, it is first of 
interest to analyse the marginally unstable states. The threshold curves shown in 
Figure 2 have symmetry: a rotation of tt radians leaves the thresholds unchanged. 
Thus, it is sufficient to show how various properties (e.g., magnetic twist and 
energy release) vary along the parts of the threshold where ai > (Figure 3). 
For ease of plotting we can convert the threshold curves to a one dimensional 
form: the filled circles and associated numbers of Figure 3 (left) represent the tic 
marks and labels for the ID threshold point axis, see Figures 7-11. Such figures 
will always be calculated using a loop of aspect ratio 20. 



Threshold Point Threshold Point 20 




Figure 4. The linear eigenfunctions, V x{x, y=0, z), for a-space points (left) and 20 
(right). The a coordinates associated with each eigenfunction are on the unstable side 
of the threshold point number (Figure 3). Cartesian coordinates are used, hence, the 
X-axis is equivalent to the radial axis. 
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Threshold Point 40 Threshold Point SO 




Figure 5. The linear eigenfunctions, Vx{x, y=0, z), for a-space points 40 (left) and 60 
(right). 

Firstly, wc plot the linear eigenfunctions for a selection of marginally unstable 
a-space points that follow the instability threshold. The location of these points 
can be determined from the threshold point number given at the top of each plot 
in Figures 4-6. The eigenfunctions of the marginally unstable modes are stongly 
radially confined; that is, there is almost no disturbance beyond the loop radius 
(i?3 = l). This contrasts sharply with the situation for loops with net current 
(Browning et al., 2008; BBV2010), in which the eigenfunction generally extends 
well into the potential envelope. 



Threshold Point 80 Threshold Point 100 




Figure 6. The linear eigenfunctions, Vx{x,y—0,z), for a-space points 80 (left) and 
100 (right). 

3.1. Instability threshold and critical twist 

Following BBV2010, we look for a single twist-related parameter that takes on 
a critical value whenever the loop reaches the threshold (Appendix B shows 
the twist profiles for a selection of loop configurations, stable and unstable). It 
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has often been postulated that instabiUty can be identified witli a single critical 
twist value irrespective of the detailed field profiles. The average twist can be 
calculated in several ways; 



C LBe{r)dr 
/J'^^ rB,Xr)dr 



Jo 
1 

-R3 JO 

1 



R3 



LBg(r) 
rB,{r) 



dr 



^Rl JO 



ZTrr -— dr 

rB,{r) 



(21) 
(22) 
(23) 



Equation 23 is the average twist weighted by area, while 21 and 22 have been 
used by VcUi, Einaudi, and Hood (1990) and Baty (2001). Note, Equation 21 
can be calculated analytically, see Appendix A. (<p)^^ denotes the average twist, 
weighted by area, over the core, outer layer and current neutralisation layer. The 
tilde {^) and hat (A) symbols are used to indicate the other equations. 



35 
30 

25 
20 
15 
10 
5 

-5 
-10 




20 



40 60 80 
Threshold Point 



100 



Figure 7. The variation in the loop average twist along the ID representation of the 
instability threshold (i/7?3 = 20). The solid lines were calculated from Equation 23; 
the long dashed from Equation 22 and the short dashed from Equation 21. The twist 
values are plotted in units of tt. 



None of the twist averages (Figure 7) is invariant around the whole threshold 
curve, although ~ Stt (Equation 22) for the majority of threshold points. 

This value is approximately in line with the oft-quoted result of 2.497r, the critical 
twist for a loop of aspect ratio 10 (Hood and Priest, 1981). Each threshold point 
has a radial twist profile; these profiles feature reversed twist until around point 
60, where the profile becomes single signed. After this point, the three average- 
twist plots converge to values between Stt and IOtt. At higher threshold points, 
the plots diverge and for Equation 22 and 23 the averages increase sharply. 

Finally, we consider the proposal of Malanushcnko et al. (2009), that a critical 
value of normalised helicity (equivalent, in our terms, to the normalised loop 
helicity, K/ip'^, over the range 0- R3) indicates instability onset. Figure 8 shows 
that the normalised helicity is certainly not the same for every threshold point, 
even if ai and a2 have the same sign. 
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a 1 




20 40 60 80 100 
Threshold Point 



Figure 8. The variation in K/tp^ (over the range O--R3) along the ID representation 
of the instability threshold (i / R-i =20). 



3.2. Path to instability 

BBV2010 used a random walk process to simulate a loop being twisted by 
turbulent photospheric motions. In other words, a loop performed a sequence 
of fixed-length steps of random direction within a-space until the instability 
threshold was crossed. We will follow this process for zero-net-current loops too, 
however, we will also employ spatially correlated random walks. This is to allow 
the correlation between the inner and outer parts of the loop to be varied. In 
particular, it is more likely that the twisting will be fairly uniform across the 
loop (i.e., the change in ai is similar to the change in 02). 

When a loop begins its random walk (i.e., when it emerges from beneath the 
photosphere) it is assigned a random starting position within the stable region 
of a-space equilibria (i.e., the loop may have some initial twist). It is more likely 
however, that the initial twist will be small and that the initial value of a2 is 
similar to (or correlated with) the initial ai value. Furthermore, the change in a- 
coordinates that occur whenever the loop steps through a-space, in response to 
photospheric driving, should also be correlated. The initial ai coordinate of the 
walk is chosen from a normal distribution centred on zero. A standard deviation 
is chosen such that the probability of the initial ai value representing an unstable 
configuration is negligible. Similarly, the initial a2 coordinate is chosen such that 
the mean is the initial ai coordinate. 

The step values, 5ai and 5a2, are determined by assuming a step size. A, 
and 5ai w 5a2- Hence, 5ai is also chosen from a normal distribution, but this 
time the mean is ^ and 5a2 is chosen such that the mean is 5ai. As the 

v2 

standard deviation of the normal distribution used to select 5ai and 5a2 is 
decreased, the range of threshold crossings narrows. In other words, the walks 
follow the ai = a2 hue more closely. A standard deviation of 0.1 will be used for 
the simulations presented in this paper, since this value restricts the threshold 
crossings to points where the twist is single-signed, see Figure 9. Correlated walks 
therefore are predisposed to maintaining the realism of loop configurations, since 
it is expected that in general photospheric motions do not create loops that have 
reversed twist. 
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Figure 9. The stability region for a loop of aspect ratio 20 is demarcated by instability 
thresholds (solid lines) and reversal lines (short dashed lines). The loop configura- 
tions along the threshold have single-signed twist (black) or reversed twist (gray). The 
relaxation line (long dashed) comprises the points within the stability region where 
ai — a2. 



Figure 9 shows that a loop might cross a 5^ reversal hne before it reaches 
the instabihty threshold. If this happens, the loop is discarded and the simu- 
lation resumes with a new loop that has a stable a-configuration. Once a loop 
reaches the instability threshold, it becomes linearly unstable. At this point, the 
field releases energy and transitions to a lower-energy state defined by Taylor 
relaxation: helicity is conserved and the a-profilc relaxes to a single value. 

3.3. Energy release calculation 

We allow each loop of an ensemble of 10^ loops to undergo a single relaxation 
(BBV2010). Initially, a loop starts from an assigned stable state. The field profile 
then undergoes a random walk (which may or may not be correlated) until it 
crosses the instability threshold; whereupon, the loop relaxes and the profile 
transitions to the relaxation line (ai =012). The relaxation a (oix) will, of course, 
vary depending on where the threshold was crossed; is found by helicity 
conservation (Taylor, 1974; Heyvaerts and Priest, 1984; Taylor, 1986; Browning 
and Van der Linden, 2003). In mathematical terms, wc find the roots of the 
following equation: 

Kia,)~Kiaa,a^2) = 0, (24) 

where an and 0^2 are the coordinates of the instability threshold crossing (con- 
servation of axial flux is assured through the normalisation -0* = 1). The helicity 
can be expressed as follows: 

K = 2Lr&ti^dr, (25) 

where /(r) is the axial current and L is the loop length (Finn and Antonsen, 
1985). Axial flux is also conserved. The full expressions for helicity and energy 
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are given in Appendix A. The energy difference between the unstable and relaxed 
states can be calculated thus: 

5W = W{aa.ar2)~W{a^). (26) 

This is the relaxation energy: the energy that is liberated from the magnetic 
field during the event. How much of this energy is converted to heat depends on 
the plasma response; thus, the relaxation energy represents the upper limit of 
the energy available for plasma heating. 

In BBV2010, a loop with net current was relaxed such that the a-profile be- 
came invariant over the range O--R4. Hence, the relaxed state always represented 
a threefold radial expansion of the threshold state (i.e., from i?3 to R4), the 
relaxation encompassed both the loop and the potential envelope. Numerical 
simulations (Browning et al., 2008) indicate that this is a good model for loops 
with net current. However, for loops with zero net current, the instability is more 
radially confined and the reconnection activity is correspondingly localised; it 
does not extend to the outer boundary (Hood, Browning, and Van dor Linden, 
2009; Bareford et al., 2011). 

We therefore consider that the relaxation radius, Rx^ can be anywhere in 
the range i?3(= 1) < Rx < Ri{^ 3). If Rx = Ra, we have complete relaxation as 
previously considered (Browning and Van der Linden, 2003; Browning et al., 
2008; BBV2010); otherwise relaxation is localised, a is constant between and 
Rx and the fields in the remaining envelope (where a = and Rx <r< R4) 
are fixed so that they do not change during relaxation; this maintains current 
neutralisation, albeit via an infinitely thin current-neutralising surface. Axial 
fiux is conserved, such that i/; (over 0-Rx) of the threshold state is equal to 
-0 (over 0-Rx) of the relaxed state. K /ip"^ is conserved in an identical manner 
(in our previous work, conservation was always over 0-i?4 and since the total 
axial flux was normalised to 1, conserving K ji^P' was identical to conserving K\ 
Likewise, the energy release is the energy of the threshold state over Q-Rx minus 
the energy of the relaxed state over the same radial range. In fact, the energy of 
the remaining potential envelope is unchanged, so that the energy release could 
also be taken over the entire volume (0-i?4); similarly, the envelope has zero 
helicity before and after relaxation. 

4. Distribution of energies and coronal heating considerations 

We now proceed to the main task of our work, which is to calculate the distri- 
bution of magnitudes of the sequence of heating events generated by random 
photospheric driving. First, we show how various properties vary along the 
instability threshold. 

4.1. Helicity and energy 

The left panel of Figure 10 plots total helicities of the threshold states. A to- 
tal helicity (or fiux) is one calculated over the range 0-i?4, i.e., the loop and 
envelope. None of the threshold states have sufficient helicity for the relaxed 
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Figure 10. Total helicity (left) and total (dimensionless) magnetic energy (right) along 
the ID representation of the instability threshold {L / R3 = 20). 



state to feature helical modes (Taylor, 1986) and so all relaxed states are cylin- 
derically symmetric. The left panel of Figure 11 confirms that each threshold 
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Figure 11. (left) and energy release (right) along the ID representation of the insta- 
bility threshold [L / R3, = 20). These properties have been calculated for two relaxation 
radii Rx = Ri (dashed) and Rx — R4 (solid). When Rx ~ R4, ctx is of 0(10^^) and so 
the corresponding plot appears very close to the = line. 

state corresponds to a relaxed state. Both of the graphs in Figure 11 feature 
two plots; the dashed line represents minimum relaxation {Rx=R3) and the 
solid line represents full relaxation [Rx = Ri). The right panel shows that, in 
general, 6W* is affected by Rx (although, there is one part of the threshold 
where the energy release is insensitive to relaxation radius); hence, the energy 
distributions that appear in the next section are calculated for minimum and 
maximum relaxation radii. 

The energies shown in Figure 10 (right) and Figure 11 (right) are given as di- 
mensionless quantities; BBV2010 derived the following expression for calculating 
a dimensional energy, 

SW = 81— RlB^ SW* , (27) 
Mo 

where Rc is the loop radius in the corona and Be is the mean axial field in 
the corona. Assuming typical values {Rc = l Mm and 5c = 0.01 T), we obtain 
dimensional energy values of 6 x 10^^ S W* J = 6 x 10^^ S W* erg. Thus, the top 



0.055 



k 0.03 
«= 0.025 
0.02 
0.015 
0.01 
0.005 




I ' ' ' ^ 

20 40 60 80 100 

Threshold Point 



SOLA: paper.tex; 18 January 2013; 17:37; p. 15 



M. R. Bareford, P. K. Browning, and R. A. M. Van der Linden 



end of the SW* scale (« 0.05) for i?^ = 3 is equivalent to 3 x 10^® erg. This is 
in the microflare range, but nanoflare energies will be obtained for weaker fields. 

4.2. Flare energy distribution 

An expression for the energy flux is derived by considering the loops in the 
ensemble as spatially separated but flaring simultaneously. All the energy in- 
put from the photosphere is dissipated, in a long-term time-average over many 
events, since the instability threshold limits the accumulation of stresses within 
the coronal magnetic field. The energy flux, F, is thus; 

10^ 

where N is the average number of steps taken to reach the threshold and r is 
the time taken to complete each step in the random walk. Similar expressions 
exist in the literature based on random photospheric twisting (Sturrock and 
Uchida, 198f; Berger, 1991; Zirker and Cleveland, 1993; Abramenko, Pevstov, 
and Romano, 2006). 

In other words, r is the time taken for a to change by A / R3', we may estimate, 
a timescale for this process as follows. Based on axial values, A corresponds to a 
change in magnetic twist Scj)^ (L/2)(A /R3); taking L/i?3 = 20 and A=: 1 gives 
6(j)~ 10. If this is caused by photospheric twisting motions of magnitude vg 
for a time interval t, we find t ^ {6(j))Rf /vg, where Rf is the footpoint radius 
(at the photosphere). With typical values of i?/ = 200km and ?;6i = lkms^^, 
we obtain t « 2000 s; note that this is consistent with the expected correlation 
time for photospheric motions; granule lifetimes are of the order 10'^ s (Zirker 
and Cleveland, 1993). Hence r has a linear relationship with the loop length, 
r=100(L [Mm]). Applying the previously used values for Be and Rc, gives a 
dimensional flux of (10^/A^r) ^ SW* erg cm~^ s~^. This result is applicable to 
Active Regions (a value for the Quiet Sun can be obtained by setting Be = 0.001 
T; this simply lowers the multiplier (10*) by 2 orders of magnitude. 

The energy distributions given below are each derived from a loop ensemble, 
that is a collection of 10^ loops flaring simultaneously. Since each loop relaxes 
only once we can sidestep the complications that come with allowing loops to 
survive many relaxations: for example, a loop may shrink or implode after flaring 
(Janse and Low, 2007), and a different instability threshold would have to be 
applied should the aspect ratio be altered as a consequence. 

4-. 2.1. Distribution of "nano flares" 

Examination of Figure 12 yields three key points. Firstly, the total energy re- 
leased increases with aspect ratio, but the average step count, N, decreases. 
This is to be expected since loop volume increases with aspect ratio, whereas 
the size of the stability region shrinks, see Figure 2. Secondly, as indicated before 
(Figure 11, right), increasing the relaxation radius increases the energy released. 
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And thirdly, correlated walks mean higher step counts, however, whether or not 
there is also an increased energy release depends on the relaxation radius. 

If Rx = R3i= 1), the energy release from correlated walks is reduced compared 
to the uncorrelated distributions; whereas, complete relaxation, Rx = Ri{=-i), 
leads to an increased energy release. This less-than-straightforward point is 
consistent with the plot that shows the variation in energy release along the 
threshold for both values of R^, see Figure 11 (right). A correlated walk would 
favour crossings around threshold point 90; when R^ ~ R3 the energy release is 
almost at its lowest for this part of the threshold, whereas the opposite is the 
case when R^ — Ra- This is also true for the thresholds applicable to loops of 
aspect ratio 10 and 30. 

For loops of aspect ratio 10, correlated walks produce distributions that have 
high-energy cut-offs - this feature is an artifact of the simple two-a model. It is 
caused by the fact that when L /R3 = 10, the relaxation line intersects the 
reversal line before the instability threshold. 
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Figure 12. Flare energy distributions for a 10^ loop ensemble, with each loop undergo- 
ing one relaxation event. The relaxation radius (Rx) associated with each event is R3 for 
the top two and R4 for the bottom two. The plots on the left correspond to uncorrelated 
random walks, those on the right to correlated driving. The distribution curves are 
colour-coded according to aspect ratio: red denotes L / R3 — 10, black i/i?3 = 20 and 
blue i/i?3 = 30. In addition, two properties are displayed for each plot, "^SW*, the 
total energy release (dimensionless) and TV, the average number of steps taken to reach 
the threshold. 



When one calculates the dimensional heat fluxes (Equation 28) one finds that 
flux is weakly dependent on aspect ratio. Further examination reveals that any 
dependence on aspect ratio can only come from ^ SW* , which is determined by 
the coordinates of the instability threshold. 5W* incorporates a length factor in 
units of the loop radius, i.e., (L/Rc) 5w* ~ dW* , where Sw* is the dimensionless 
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energy release per unit of dimensionless length. Substituting the full expression 
for the step time (r = {\/2){L/ve){Rf /Rc), Section 4.2) into Equation 28 gives 

10^ 

again, applying previously used values, this simplifies to = (10^/7V) ^ Sw* 
erg cm~^ s~^. The length terms cancel and the ratio Rc/Rf is effectively a 
constant. 

For distributions derived from uncorrelated walks and minimal relaxation 
(Rx — R-s), f w3-4xl0^ erg cm~^ s~^. Using correlated walks instead, dimin- 
ishes the fluxes to 0.9-2x10^ erg cm~^ s~^. Increasing the relaxation radius 
to i?4 will reverse this reduction and yield Fpa7-I0xl0^ erg cm~^ s'-^. This 
last result is also true for distributions based on uncorrelated walks and full 
relaxations. When R^ = R4 correlated walks do lead to higher energy releases, 
however, these walks are longer and have higher step counts, which means the 
flux remains roughly constant. 

Finally, in Figure 13 we show the logs of the flare energy distributions pre- 
sented in Figure 12. The distributions calculated from uncorrelated walks give 
log plots that almost match the critical gradient for coronal heating. Although 
the log plots of the correlated (Gaussian-shaped) distributions do not follow 
power-laws, we include these results for completeness. 
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Figure 13. The logarithm of the flare energy distributions for a IC loop ensemble, with 
each loop undergoing one relaxation event. The presentation of these plots follows the 
scheme used for Figure 12. The grey diagonal line in each plot is there for comparison; 
it has a gradient equal to the critical gradient for coronal heating, -2. 
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5. Summary and conclusions 

We have investigated tlie distribution of energy releases in an ensemble of coronal 
loops driven by random photosphcric footpoint motions, using Taylor relaxation 
theory. The twisting has been assumed to be localised within the loop cross 
section, so that the loop is always without net current (the azimuthal field 
vanishes at - and beyond - the loop boundary). This means we are genuinely 
studying individual loops, rather than (unrealistically) allowing the potential 
envelope outside the loop to be twisted as the loop evolves, as in previous work 
(BBV2010). 

A relaxation event is triggered whenever the loop becomes unstable to an 
ideal kink instability; during the nonlinear phase, current sheets form and subse- 
quent rapid reconnection occurs. A distribution of events is built up by allowing 
loop equilibria to evolve through a random walk, representing the effects of 
turbulent photosphcric footpoint motions, until the linear stability threshold is 
reached. The effectiveness of the relaxation approach is that energy release is 
easily calculated for a wide family of profiles, this is extremely difficult with 3D 
numerical simulations. Furthermore, relaxation theory becomes a better repre- 
sentation for very high conductivities, which cannot be accessed by present day 
simulations. This approach can, of course, be extended to more complex field 
models than the simple cylindrical coronal loop models used here. The energy 
fluxes obtained are sufficient for coronal heating; the fluxes agree with the oft- 
quoted estimates of Withbroe and Noyes (1977). The heating flux is larger for 
photosphcric motions with higher temporal correlation; within our model, this 
is represented by larger step sizes for the random walk (so that the loop is 
coherently twisted until it reaches instability, rather than randomly twisting 
and untwisting many times). We thus show that dissipation within loops could 
be sufficient for coronal heating, but in reality, this is likely to form only part of 
the coronal heating input. Topological complexity arising from, for example, the 
discrete nature of photosphcric flux sources (Priest, Heyvaerts, and Title, 2002) 
and braiding motions (Parker, 1972), will also play a strong role. 

A distribution of heating events, or nanoflares, is obtained, for a variety of 
conditions. For the case of spatially uncorrelated twisting motions, in which the 
motions may vary strongly across the loop cross section, a power law distribution 
of energy versus occurrence frequency is obtained, with a slope slightly steeper 
than the critical value of -2 required for nanoflare heating to be effective (Hudson, 
1991). For strongly correlated twist motions, in which the twist in the outer part 
of the loop is close to that in the inner core, a peaked energy distribution is 
obtained, with almost Gaussian shape. The former case reflects the distribution 
of available energies around the instability threshold, whereas as the latter is 
mainly determined by the allowable range of twist profiles. It should be noted 
that these distributions (Figures 12-13) are obtained for an ensemble of identical 
loops: in reality, much broader distributions will result due to variations in axial 
field strengths and photosphcric driving. The true nanoflare distribution is a 
convolution over more than one parameter. 

The effect of loop aspect ratio has been considered and been found to have 
little impact on energy flux. The higher volume of large aspect ratio loops is 
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counteracted by the smaller stability region (instability occurs at lower a values). 
As the aspect ratio is increased beyond 30, we expect the stability region to 
reduce by smaller and smaller amounts. In other words, the region will converge 
to a minimum area. This has been shown for constant-a loops, see Figure 4 
of Browning and Van der Linden (2003). Hence, assuming this expectation is 
verified, the energy flux will be independent of aspect ratio (Equation 29), 
assuming that the same axial field strength is applied to all members of the 
loop ensemble. Presumably, there is a dependence between loop size and \Bc\, so 
an ensemble that features some distribution of field strengths will still depend 
(albeit indirectly) on the aspect ratio. 

Contrary to the Bz profiles of Appendix B, the axial field at the loop foot- 
points should not change during the random walk or during relaxation. The 
reason for this discrepancy is that preservation of the footpoint axial field intro- 
duces a z dependency - the field becomes two dimensional. However, if the length 
of the loop exceeds its radius, a ID field approximation - such as the one used 
by the model presented here - still remains adequate for a substantial portion of 
the loop. Zweibel and Boozer (1985) and Browning and Hood (1989) show that 
the z dependence is confined to thin boundary layers near the footpoints. Hence, 
the difference in energies for loops represented by ID and 2D fields is negligible 
especially ii L / R3 > 10 (see also Lothian and Browning, 2000; Robertson, Hood, 
and Lothian, 1992). Dalmasse, Browning and Bareford (2011) have investigated 
a simpler loop, having just a core and outer layer (with a conducting wall at 
r = 1 ) , by calculating the energy releases according to Taylor relaxation for a 
representative sample of threshold configurations. This was done using both ID 
and 2D fields, with the latter maintaining the axial field at the footpoints. The 
resulting energy releases differ by less than 1% between the ID and 2D cases. 

The results presented here are based on a loop model that has a thin current- 
neutralising layer (this approximates to a current sheet), in which the fields 
discontinuously change at the loop edge. The main reason for this choice is so 
that the fields inside the loop are close to the previously-studied two-a model 
(BBV2010), and thus a comparison can be made with previous work. Also, such 
fields correspond to twisting within an isolated fiux source, whilst the fiux which 
surrounds the loop in the corona originates from untwisted separated sources. 
Interestingly, the ideal instability threshold in this case is very similar to that 
found for a close-fitting conducting wall at the loop edge, as originally used by 
Browning and Van der Linden (2003). This is because the thin current layer 
forces unstable perturbations to vanish (almost) at the loop edge. In numerical 
simulations, the choice of a thin current layer has consequences in allowing 
resistive modes to be significant; although for realistic values of the resistiv- 
ity (unattainable in simulations) the growth rate of such modes is extremely 
slow. Preliminary studies have also been undertaken with a thicker current- 
neutralising layer. In this case, a closed stability threshold curve can be obtained, 
and the results are more similar to previous work (BBV2010). 

One important consequence of considering loops with zero net current is that 
the reconnection activity tends to be localised near the loop and thus relaxation 
is likely to be incomplete (rather than including a large part of the surrounding 
potential field). We consider here two limiting cases: localised relaxation, in 
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which only the loop volume relaxes to a minimum energy (constant-a) state, and 
the surrounding potential envelope remains unaffected; and complete relaxation, 
in which the loop and the potential envelope relax out to the external boundary. 
The latter is clearly the true minimum energy state. Numerical simulations 
(Barcford et al., 2011) indicate an intermediate situation, but somewhat closer 
to the completely localised relaxation. In fact, the loop reconnects with some of 
the surrounding axial field, but only to a limited extent. This is an important 
issue for understanding relaxation in the Sun, where the extent of relaxation is 
not defined by conducting walls: in contrast with laboratory plasmas (Taylor, 
1974). This is discussed more fully in a companion paper (Barcford et al., 2011). 
In general, complete relaxation naturally gives larger energy releases, but the 
choice of relaxation radius does not strongly affect the distribution of heating 
events. Future work will use numerical simulations to explore the transition to 
instability, and the effects of continual driving. 
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Appendix 

A. Expressions for loop properties 

Expressions for some key quantities {{if), K and W) are given here. For compactness, 
tliese are given only for ai 7^ and 02 ^ 0, while special cases (e.g., a\ — 0) must be dealt 
with separately. Expressions for constant-o? fields can be recovered by setting ai = 02, 
which gives more familiar formulae. The superscripts and subscripts that accompany 
each quantity term denote the upper and lower radial bounds over which the quantity 
is calculated. The vacuum permeability, fio, used in the magnetic energy expressions, 
is set to 1. 

A.l. Average magnetic twist 
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A. 2. Magnetic helicity 
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A. 3. Magnetic energy 
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B. Magnetic field profiles for a selection of a-space points 

The magnetic axial field, Bz{r), azimuthal field, Be{r), and magnetic twist, 0(r), 
profiles are presented for a selection of stable and unstable loop configurations, see 
Fig. 14. 
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Figure 14. The stability region for a loop of aspect ratio 20 is demarcated by instability 
thresholds (solid lines) and reversal lines (short dashed lines). The relaxation line 
(long dashed) comprises the points within the stability region where ai=a2- Stable 
configurations are indicated by empty circles and unstable ones by filled circles. 
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a, = 0.8, ttj = 2-0 (stable) 



a, = 0.5, = -3.0 (stable) 



0.08 
0.06 
0.04 
0.02 


-0.02 



0.08 I 

0.06 - 

0.04 

0.02 - 

- - 

-0.02 - 



0.5 1 1.5 2 2.5 3 



0.5 1 1.5 2 2.5 3 



a, = -1.0, = 2.0 (stable) 



0.08 I 

0.06 

0.04 -- 

0.02 

■ -. 
-0.02 



= 3.4, Kg = 0.8 (unstable) 



0.08 
0.06 
0.04 
0.02 


-0.02 



0.5 1 1.5 2 2.5 3 



0.5 1 1.5 2 2.5 3 



Figure 15. The B~ (solid) and Be (dashed) profiles for some of the configurations (3 
stable, 1 unstable) located on Figure 14. 
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Figure 16. The magnetic twist (in units of tt) profiles for some of the configurations 
(3 stable, 1 unstable) located on Figure 14. 



The empty circle on the a.\ =a2 line in Fig. 14 is the relaxed state of the unstable 
configuration identified by the filled circle on the threshold. The radius of the relaxed 
state is 1.5. 
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Figure 17. The B- (solid) and Be (dashed) profiles for the threshold and relaxed 
configurations located on Figure 14. 
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Figure 18. The magnetic twist (in units of tt) profiles for the threshold and relaxed 
configurations located on Figure 14. 
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